
###########################################################
##### Haiti elite network project  		          			#####
##### Oster (2019) sensitivity analysis            		#####
##### 2021 mar 03                   									#####
###########################################################


library(robomit)


#####
## read in data
#####

fam <- read.csv('01_Data/02_Clean/fam.csv')
all <- read.csv('01_Data/02_Clean/allfams.csv')



#####
## estimate R-square from original estimations
#####

all$clust <-  as.factor(all$fastgreedy)
fam$clust <-  as.factor(fam$fastgreedy)

## biz elite

mod8 <- lm(coup ~ bonw_02_wnind_st + value_log_st + cshare2_st + noconsump +
             ref_lib_st + pci_value_st + divis_ln_st + bulk_ln_st + time_con_st + 
             pol + mil + syrian + immig + nind_log, data = fam)
mod8$se <- coeftest(mod8, vcov = vcovHC, type = "HC1")[,2]
mod8$se_clust <- coeftest(mod8, vcov = cluster.vcov(mod8, fam$clust))[,2]
inc.obs <- complete.cases(fam[,colnames(mod8$model)])
mod8$nclust <- length(unique(fam$clust[inc.obs==T]))

mod8b <- lm(coup ~ bonw_02_wnind_st + value_log_st + cshare2_st + noconsump +
              ref_lib_st + pci_value_st + divis_ln_st + bulk_ln_st + time_con_st + 
              pol + mil + syrian + immig + nind_log + clust, data = fam)
mod8b$se <- coeftest(mod8b, vcov = vcovHC, type = "HC1")[,2]
mod8b$se_clust <- coeftest(mod8b, vcov = cluster.vcov(mod8b, fam$clust))[,2]
inc.obs <- complete.cases(fam[,colnames(mod8b$model)])
mod8b$nclust <- length(unique(fam$clust[inc.obs==T]))

## all elite

mod3 <- lm(coup ~ bonw_02_wnind_st + biz + pol + mil + syrian + immig + nind_log, data = all)
mod3$se <- coeftest(mod3, vcov = vcovHC, type = "HC1")[,2]
mod3$se_clust <- coeftest(mod3, vcov = cluster.vcov(mod3, all$clust))[,2]
inc.obs <- complete.cases(all[,colnames(mod3$model)])
mod3$nclust <- length(unique(all$clust[inc.obs==T]))

mod3b <- lm(coup ~ bonw_02_wnind_st + biz + pol + mil + syrian + immig + nind_log + clust, data = all)
mod3b$se <- coeftest(mod3b, vcov = vcovHC, type = "HC1")[,2]
mod3b$se_clust <- coeftest(mod3b, vcov = cluster.vcov(mod3b, all$clust))[,2]
inc.obs <- complete.cases(all[,colnames(mod3b$model)])
mod3b$nclust <- length(unique(all$clust[inc.obs==T]))



## importers sample  

  ## plot betas over range of delta [0.1,1]

tab <- data.frame('beta_stars' = rep(NA, 10),
                  'lo95' = rep(NA, 10),
                  'up95' = rep(NA, 10))

d <- seq(0.1, 1, 0.1)

for (i in 1:10){
  delta <- d[i]
  t <- o_beta(y = "coup", # dependent variable
         x = "bonw_02_wnind_st", # independent treatment variable
         con = "value_log_st + cshare2_st + noconsump +
              ref_lib_st + pci_value_st + divis_ln_st + bulk_ln_st + time_con_st + 
              pol + mil + syrian + immig + nind_log + clust", # other control variables
         delta = delta, # delta
         R2max = 1.3*summary(mod8b)$r.squared, # maximum R-square
         type = "lm", # model type
         data = fam) # dataset
  tab$beta_stars[i] <- t[1,2]
  
  t <- o_beta_boot(y = "coup", # dependent variable
                   x = "bonw_02_wnind_st", # independent treatment variable
                   con = "value_log_st + cshare2_st + noconsump +
              ref_lib_st + pci_value_st + divis_ln_st + 
              bulk_ln_st +
              time_con_st +
              pol + mil + syrian + immig + nind_log + clust", # other control variables
                   R2max = 1.3*summary(mod8b)$r.squared, # maximum R-square
                   delta = delta, # delta
                   sim = 1000, # number of simulations
                   obs = 200, # draws per simulation
                   rep = T, # bootstrapping with or without replacement
                   type = "lm", # model type
                   useed = 123, # seed
                   data = fam) # dataset
  
  tab$lo95[i] <- unlist(sort(t['beta*']))[25]
  tab$up95[i] <- unlist(sort(t['beta*']))[975]
}
tab$delta <- d
tab$beta_stars <- unlist(tab$beta_stars)

pdf("03_Figures/sens_oster_fam.pdf")
ggplot(tab, aes(x = delta, y = beta_stars)) +
  geom_point(aes(x = delta, y = beta_stars)) +
  theme_minimal() + 
  geom_errorbar(aes(ymin = lo95, ymax = up95)) +
  ylim(-0.4,0.4) +
  xlab("Adjusted Coefficient")
dev.off()


o_delta(y = "coup", # dependent variable
        x = "bonw_02_wnind_st", # independent treatment variable
        con = "value_log_st + cshare2_st + noconsump +
              ref_lib_st + pci_value_st + divis_ln_st + 
              bulk_ln_st +
              time_con_st +
              pol + mil + syrian + immig + nind_log + clust", # other control variables
        R2max = 1.3*summary(mod8b)$r.squared, # maximum R-square
        beta = 0, # beta
        type = "lm", # model type
        data = fam) # dataset



## all elite sample  

  ## plot betas over range of delta [0.1,1]

tab <- data.frame('beta_stars' = rep(NA, 10),
                  'lo95' = rep(NA, 10),
                  'up95' = rep(NA, 10))

d <- seq(0.1, 1, 0.1)

for (i in 1:10){
  delta <- d[i]
  t <- o_beta(y = "coup", # dependent variable
              x = "bonw_02_wnind_st", # independent treatment variable
              con = "biz + pol + mil + syrian + immig + nind_log + clust", # other control variables
              delta = delta, # delta
              R2max = 1.3*summary(mod3b)$r.squared, # maximum R-square
              type = "lm", # model type
              data = all) # dataset
  tab$beta_stars[i] <- t[1,2]
  
  t <- o_beta_boot(y = "coup", # dependent variable
                   x = "bonw_02_wnind_st", # independent treatment variable
                   con = "biz + pol + mil + syrian + immig + nind_log + clust", # other control variables
                   R2max = 1.3*summary(mod3b)$r.squared, # maximum R-square
                   delta = delta, # delta
                   sim = 1000, # number of simulations
                   obs = 200, # draws per simulation
                   rep = T, # bootstrapping with or without replacement
                   type = "lm", # model type
                   useed = 123, # seed
                   data = all) # dataset
  
  tab$lo95[i] <- unlist(sort(t['beta*']))[25]
  tab$up95[i] <- unlist(sort(t['beta*']))[975]
}

tab$delta <- d
tab$beta_stars <- unlist(tab$beta_stars)

pdf("03_Figures/sens_oster_all.pdf")
ggplot(tab, aes(x = delta, y = beta_stars)) +
  geom_point(aes(x = delta, y = beta_stars)) +
  theme_minimal() + 
  geom_errorbar(aes(ymin = lo95, ymax = up95)) +
  ylim(-1.1,3) +
  xlab("Adjusted Coefficient")
dev.off()

  ## estimate delta needed to bring coefficient to 0

o_delta(y = "coup", # dependent variable
        x = "bonw_02_wnind_st", # independent treatment variable
        con = "biz + pol + mil + syrian + immig + nind_log + clust", # other control variables
        R2max = 1.3*summary(mod3b)$r.squared, # maximum R-square
        beta = 0, # beta
        type = "lm", # model type
        data = all) # dataset



